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Brine Treated Biomass and the Adsorption of Dyes 
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Department of Industrial Management and Technology, University of Piraeus, Piraeus 18534, Greece 


Abstract: Brine is a solution of salt (usually sodium chloride) in water. In different contexts, brine may refer to salt solutions ranging 


from about 3.5% (a typical concentration of seawater, or the lower end of solutions used for brining foods) up to about 26% (a typical 


saturated solution, depending on temperature). Adsorption onto activated carbon is the most widespread technology for the removal 


of pollutants from water and wastewaters. In this study, continuous fixed-bed-column systems were investigated. The adsorbents 


which authors use are: spruce (Picea abies) untreated, spruce modified by autohydrolysis. The column systems were filed with 


biomass at various initial dye concentrations, flow rates and bed-depths. The column kinetics of MB (Methylene Blue) adsorption on 


spruce (Picea abies) untreated, spruce modified by autohydrolysis was simulated. Economies arise when the facility that can use such 


adsorption materials is near a source of a lignocellulosic waste as agricultural residues, thus saving transportation cost and 


contributing to industrial ecology at local level. 


Key words: Adsorption, desorption, column studies, pretreatment, brine. 


1. Introduction 


Brine is a solution of salt (usually sodium chloride) 
in water. In different contexts, brine may refer to salt 
solutions ranging from about 3.5% (a typical 
concentration of seawater, or the lower end of 
solutions used for brining foods) up to about 26% (a 
typical saturated solution, depending on temperature). 
Brine is a liquid, and 0 °F was initially set as the zero 
point in the Fahrenheit temperature scale, as it was the 
coldest temperature that Daniel Gabriel Fahrenheit 
could reliably reproduce—by freezing brine. At 
100 °C (373.65 K, 212 °F), saturated sodium chloride 
brine is about 28% salt by weight i.e. 39.12 g salt 
dissolves in 100 mL of water at 100 °C. At 0 °C 
(273.15 K, 32 °F), brine can only hold about 26% salt. 
Brine is a common fluid used in large refrigeration 
installations for the transport of thermal energy from 


place to place. It is used because the addition of salt to 
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water lowers the freezing temperature of the solution 
and the heat transport efficiency can be greatly 
enhanced for the comparatively low cost of the 
material. The lowest freezing point obtainable for 
NaCl brine is -21.1 °C (-6.0 °F) at 23.3wt% NaCl. 
This is called the eutectic point [1, 2]. According to 
many review papers [3-8], low-cost adsorbents offer a 
lot of promising benefits for commercial purposes in 
the future. They could be used in place of commercial 
activated carbon for the removal of dyes in aqueous 
solutions. 

Adsorption onto activated carbon is the most 
widespread technology for the removal of pollutants 
from water and wastewaters. The disadvantage of 
activated carbon is its high cost [11]. Hence, it is of 
pivotal importance thence of low-cost substitute 
absorbents to replace activated carbons. Various types 
of untreated biomass have been reported to have a use 
in dye removal: sawdust [9, 11], wheat straw [12], 
cedar sawdust [13], rubberwood sawdust [14], kudzu 


[15], banana and orange peels [16] and palm kernel 
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fiber [11, 17], peanut husk [4], palm kernel fibre [11], 
turbinate alga [12], graphene [13], 
defatted jojoba [14] and sugar bet pulp [15]. 

Further, 


materials are used to remove dyes in water and 


Turbinaria 


numerous pretreated lignocellulosic 
wastewater: Pyrolyzed date pits [16], date stones [17] 
and Turbinaria turbinate alga. 

In this study, continuous fixed-bed-column systems 
were investigated. The adsorbents which authors use 
are: spruce (Picea abies) untreated and spruce 
modified by brine treatment. The column systems 
were filed with biomass at various initial dye 
concentrations, flow rates and bed-depths. The column 
kinetics of MB (Methylene Blue) adsorption on spruce 
(Picea abies) untreated, and spruce modified by brine 
treatment, was simulated herein, using biomass as 
control, in order to facilitate its potential use as a low 
cost adsorbent for wastewater dye removal. 
Economies arise when the facility that can use such 
adsorption materials is near a source of a 
lignocellulosic waste as agricultural residues, thus 
contributing to 


saving transportation cost and 


industrial ecology at local level. 


2. Materials and Methods 


The brine treatment process was performed in a 
3.75-L batch reactor PARR 4843. The isothermal 
autohydrolysis time was tai = 0, 10, 20, 30, 40 and 50 
min (not including the non-isothermal preheating and 
the cooling time-periods); the reaction was catalyzed 
by the organic acids produced by the pine sawdust 
itself during autohydrolysis at a liquid-to-solid ratio of 
10:1; the liquid phase volume (water) was 2,000 mL 
and the solid material dose (pine sawdust) was 200 g; 
150 rpm. The 
temperature values were T = 160 °C, 200 °C and 
240 °C, reached after t = 42, 62 and 80 min preheating 
respectively. The 


stirring speed reaction ending 


time values, autohydrolysis 
product was filtered using a Buchner filter with 
Munktell paper sheet (grade 34/N) to separate the 


liquid phase and from the solid phase. The solid 


residue was washed with water until neutral pH (the 
initial filtrate pH was 2.90-4.76 depending on the 
autohydrolysis severity). The solid residue was dried 
at 110 °C for 10 days at room temperature to reach the 
humidity of the untreated material. Then it was used 
as adsorbent. 

Continue-flow experiments were carried out on 
stainless steel columns with dimensions 15 x 2.5 cm. 
The bed height was x = 15 cm, respectively. 
The adsorbent weight was m = 32 g and 54 g, 
respectively. The pH was 8.0. The flow rates were 
fixed at approximately 20 mL-min' 
preparative HPLC 
Chromatography) pump, LaPrep P110-VWR-VWR 


International. The initial concentrations of MB were 


using a 


(High Performance Liquid 


165 mg L". To determine the concentration of MB in 
the effluent, samples of outflow were peaked at 100 


mL intervals. 


3. Results and Discussion 


A widely used continuous fixed-bed-column model 
was established by Bohart and Adams [17], who 
assumed that the rate of adsorption is controlled by the 
surface binding (through chemical reaction or physical 
interaction) between adsorbate and unused capacity of 
the solid, i.e., adsorption rate = KC C,„, where K is the 
adsorption rate coefficient, C is the adsorbate 
concentration at the solid phase at distance x, and C, is 
the unused surface adsorptive capacity at time t, 
expressed as mass per volume of bed. The material 
balance for adsorbate is given by the partial 
differential equation: 


Cj 
a ey are ba | 
z u (1) 
while the corresponding partial differential equation 
for the C, decrease is 
oC K 
ECG 2) 
Ox u 
where u is the superficial liquid velocity. These 


equations are obtained neglecting diffusion and 
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accumulation terms, assumptions that are valid in 
chemical engineering practice, provided that strict 
scale up specifications are kept in the design stage and 
successful operation conditions are kept in the 
industrial operation stage. 

The differential equations can be integrated over the 


total length x of the bed to give: 
C NG 
n| £ — ) = Inex £ Yx) —I|-K-C,-t 8) 
u 


where N (mg L’) is the initial or total adsorption 





capacity coefficient, also quoted as C„o [17]; C = 
effluent concentration (mg: L’); C; = influent 
concentration (mg: LÙ; K = adsorption rate coefficient 
(L-mg'-min');, x = bed depth (cm); u = linear 
velocity (cm-min’); and ¢ = time (min). Since 
exp(K.N.x/u) is usually much greater than unity, 


this equation can be simplified to: 


nf S21) 282K. € ø 


u 


which is commonly used by researchers, because of its 
convenience in estimating the values of parameters K 
and N through linear regression either of In[(Co/C; ) — 


1] vs. t or t vs. x when the following rearrangement is 


N-x ji C. 
C u Keel C 


adopted: 








In this rearrangement, t is the time to breakthrough, 
1.e., the time period required for concentration to reach 
a predetermined value. For using the last expression as 
a linear regression model, wastewater is passed 
through beds of varying depths, keeping constant C; 
and u, preferably at values similar to those expected to 
prevail under real conditions at full scale. 
Alternatively, it can be performed by the aid of at least 
three columns arranged in series. In such a case, 
sampling takes place at the bottom of each column 
and is measured for adsorbate concentration, making 


more frequent measurements when approaching the 


breakthrough concentration C. Finally, the time at 
which the effluent reaches this concentration is used 
as the dependent variable while x plans the role of the 
independent one. Evidently, the use of such a 
regression model implies the additional error of 
measuring the independent variable with less precision 
in comparison with the dependent. The common error 
in both models comes from the estimation of 
concentration from measuring adsorbance although 
the reference relation/curve has been structured/drawn 
in the 


concentrations the corresponding adsorbances have 


inverse mode, 1.e., for predetermined 
been measured. 

In the present work, the model of Eq. (17) has been 
used for parameter values estimation through linear 
regression to obtain numerical results comparable with 
corresponding data found for other fixed bed 
adsorption studies in literature. The non-linear form of 
this model is: 

C 


Ze 6 
l+4e” ©) 


where An-ek Nxiu; r=K-C;. 
On the other hand, Clark [18] has advanced the 
Bohart and Adams model [19] by incorporating the 


parameter n of the Freundlich adsorption isotherm: 


1 
C RE 
© "| Te de a 


where n = inverse of the slope of the Freundlich 
isotherm [12]. Finally, the Bohart and Adams model 
can be reduced for n = 2 from Clark model [19]. 

The kinetic equation used for desorption is Eq. (21): 


Cage (21) 


where C’) is the initial MB concentration of 
desorption effluent, k’ is desorption rate constant 
assuming first order desorption kinetics and £’ is 
desorption time. 


The use of salt brine—a mixture of 23 percent salt 
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and 77 percent water (the carrying capacity of slat in a 
water solution)—for melting snow and ice was 
developed in Europe. The mixture is sprayed directly 
on the road as an anti-icing agent. Although salt brine 
has the same melting characteristics of regular salt, 
because it is in brine form, it has the advantage of 
working immediately and does not have the problem 
of “bouncing and scattering” off the road because of 
its consistency, according to a study of salt brine use 
by the Vermont Agency of Transportation, Materials 
and Research Section. The agency found that brine 
was cost effective to produce, at about $0.10 per 
gallon. Additionally, other deicing chemicals can be 
mixed into the salt brine to lower the effective melting 
temperature. During the study (Table 1) in which 


experimental test sections of Vermont interstates were 


treated with salt brine and salt brine mixtures, the 
agency found that it saved an average of 24 percent of 
material usage. 

Although the data sets were not as extensive as the 
research team initially thought, the research project 
was able to produce salt savings of approximately 30 
percent and materials cost savings of approximately 
30 percent and materials cost savings of 
approximately 24 percent during the first year of 
experimentation, according to the Vermont DOT 
report. If such a cost saving could be projected 
statewide, the potential savings would be almost $1 
million annually. The research team agreed that there 
is potential for more salt and cost savings as staff 
becomes more experienced with the technology (Figs. 


1-7). 


Table 1 Fixed bed column systems for spruce modified by brine treatment. 


Concen 


sacd X Te TCC) ¢t(min) Ci O (mL/min) x (cm) m (g) N K R qo (mg/g) 
l 180 50 153 18 15 20 4948 0.000353 -0.9481 18.21 
2 180 50 145 21 15 20 6975 0.000242 -0.9683 25.67 
4 140 0 210 21 15 12.11 8461 0.000153 -0.9616 51.42 
4 160 0 156 21 15 11.7 4304 0.000694 -0.9946 27.07 
4 160 50 163 21 15 20 5834 0.000300 -0.9662 21.47 
4 180 0 162 21 15 11.6 6054 0.000257 -0.9744 38.41 
4 180 50 151 21 15 20 5666 0.000444 -0.9479 20.85 
4 200 0 155 22 15 12.3 4682 0.000866 -0.9667 28.01 
4 200 50 162 21 15 14.6 9259 0.000132 -0.9390 46.67 
4 200 50 148 21 15 20 7713 0.000157 -0.9232 28.38 
4 220 50 157 21 15 12 5840 0.000235 -0.9863 35.82 
4 200 50 160 21 15 20 9233 0.000181 -0.9466 33.97 
4 240 50 162 19 15 13 1358 0.000899 -0.9237 7.69 
4 240 50 153 21 15 12 451 0.000369 -0.7918 2.76 
4 240 50 153 21 15 12 451 0.000369 -0.7918 2.76 
8 180 0 169 21 15 12 4715 0.000212 -0.9565 28.92 
8 180 50 159 17 15 173 7089 0.000351 -0.9211 30.16 
8 180 50 141 21 15 20 5695 0.000259 -0.9799 20.96 
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Fig. 1 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg'L“', O=20 mL-min™ (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated <1 time at 180 °C for 50 min. 
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Fig. 2 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg-L", O=20 mL-min™ (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated x2 times at 180 °C for 50 min. 
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Fig. 3 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg-L", O=20 mL-min™ (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated x2 times at 160 °C for 10 min. 
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Fig. 4 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg-L", O=20 mL-min” (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated x2 times at 160 °C for 50 min. 
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Fig. 5 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg-L", O=20 mL-min” (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated <4 times at 180 °C for 0 min. 
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Fig. 6 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 
presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg-L", O=20 mL-min’ (the theoretical 
curves are according to the Bohart and Adams model); modified by brine treatment; concentrated x4 times at 200 °C for 50 min. 
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Fig.7 Column experimental data and theoretical curves of MB adsorption on modified spruce; the effluent concentration is 


presented vs. (a) the effluent volume and (b) the adsorption time; x = 15 cm, C; = 160 mg L’, O=20 mL-min™ (the theoretical 


curves are according to the Bohart and Adams model); modified by brine treatment; concentrated x8 times at 180 °C for 50 min. 


References 


[1] 


Rafatullah, M., Sulaiman, O., Hashim, R., and Ahmad, A. 
2010. “Adsorption of Methylene Blue on Low-Cost 
Adsorbents: A Review.” J. Hazard. Mater. 177: 70-80. 
El-Sayed, A. M., Mitchell, V., Manning, L. A., Cole, L., 
and Suckling, D. M. 2011. “New Sex Pheromone Blend 
for the Lightbrown Apple Moth Epiphyas postvittana.” J. 
Chem. Ecol. 37: 640-6. 

Altenor, S., Ncibi, M. C., Emmanuel, E., and Gaspard, S. 
2012. 
Properties and Adsorption Efficiencies of Caribbean Alga 


“Textural Characteristics, |§ Physiochemical 


Turbinaria turbinata and Its Derived Carbonaceous 
Materials for Water Treatment Application.” Biochem. 
Eng. J. 67 (2): 35-44. 

Liu, T., Li, Y., Du, Q., Sun, J., Jiao, Y., Yang, G., et al. 
2012. “Adsorption of Methylene Blue from Aqueous 
Solution by Grapheme.” Colloids and Surf. B: 
Biointerfaces. 90: 197-203. 

Al-Anber, Z. A., Al-Anber, M. A., Matouq, M., Al-Ayed, 
O., and Omari, N. M. 2011. “Defatted Jojoba for the 
Removal of Methylene Blue from Aqueous Solution: 
Thermodynamic and Kinetic Studies.” Desalination 276: 
169-74. 

Malekbala, M., Hosseini, S., Yazdi, S. K., and Masoudi 


[11] 


Soltani, S. 2012. “The Study of the Potential Capability 
of Sugar Beet Pulp on the Removal Efficiency of Two 
Cationic Dyes.” Chem. Eng. Res. Des. 90: 704-12. 
Theydan, S. K., and Ahmed, M. J. 2012. “Adsorption of 
Methylene Blue onto Biomass-Based Activated Carbon 
by FeCl, 
Thermodynamic Studies.” J. Anal. Appl. Pyrol. 97: 
116-22. 

Ahmed, M. J., and Dhedanb, S. K. 2012. “Equilibrium 
Isotherms and Kinetics Modeling of Methylene Blue 


Activation: Equilibrium, Kinetics, and 


Adsorption on Agricultural Wastes-Based Activated 
Carbons.” Fluid Phase Equilibr. 317: 9-14. 

Oliveira, W. E., Franca, A. S., Oliveira, L. S., and Rocha, 
S. D. 2008. “Untreated Coffee Husks as Biosorbents for 
the Removal of Heavy Metals from Aqueous Solutions.” 
J. Hazard. Mater. 152: 1073-81. 

Reffas, A., Bernardet, V., David, B., Reinert, L., 
Bencheikh Lehocine, M., Dubois, M., et al. 2010. 
“Duclaux, Carbons Prepared from Coffee Grounds by 
H3PO, Activation: Characterization and Adsorption of 
Methylene Blue and Nylosan Red N-2RBL.” J. Hazard. 
Mater. 175: 779-88. 

Baek, M. H., Ijagbemi, C. O., Se-Jin, O., and Kim, D. S. 
2010. “Removal of Malachite Green from Aqueous 
Solution Using Degreased Coffee Bean.” J. Hazard. 


[12] 


[14] 


[15] 


Execution of Measurements for Determining the Parameters Affecting the Thermochemical 91 
Treatment of Brine Treated Biomass and the Adsorption of Dyes 


Mater. 176: 820-8. 

Kyzas, G. Z., Lazaridis, N. K., and Mitropoulos, A. C. 
2012. “Removal of Dyes from Aqueous Solutions with 
Untreated Biomass as Potential Low-Cost Adsorbents: 
Equilibrium, Reuse and Thermodynamic Approach.” 
Chem. Eng. J. 189-190: 148-59. 

Freundlich, H. M. F. 1906. “Uber die adsorption in 
lösungen, Zeitschrift für Physikalische.” Chemie. 57: 
385-471. 

Langmuir, I. 1916. “The Constitution and Fundamental 
Properties of Solids and Liquids.” J. Am. Chem. Soc. 38: 
2221-95. 

Sips, R. 1948. “Structure of a Catalyst Surface.” J. Chem. 


[16] 


[17] 


Phys. 16: 490-5. 
Lagergren, S. 1898. 


geloster 


“Zur theorie der sogenannten 
stoffe, 

Vetenskapsakademiens.” Handlingar 24: 1-39. 

Ho, Y. S., Ng, J. C. Y., and McKay, G. 2000. “Kinetics 

of Pollutants Sorption by Biosorbents: Review.” Sep. 

Purif. Methods 29: 189-232. 

Weber, W. J., and Morris, J. C. 1963. “Kinetics of 

Adsorption on Carbon from Solution.” J. Sanit. Eng. Div. 

Am. Soc. Civ. Eng. 89: 31-60. 

Bohart, G., and Adams, E. N. 1920. “Some Aspects of the 

Behavior of Charcoal with Respect to Chlorine.” J. Am. 

Chem. Soc. 42: 523-44. 


adsorption Kungliga Svenska 


Journal of Environmental Science and Engineering B 7 (2018) 92-102 
doi:10.17265/2162-5263/2018.03.002 





PUBLISHING 


Changing of Properties of Unsaturated Compacted 


Bentonite due to Hydration Effort 


Tomoyoshi Nishimura’ and Junichi Koseki’ 
1. Department of Civil Engineering, Ashikaga University, Ashikaga, Tochigi 326 8558, Japan 
2. Department of Civil Engineering, The University of Tokyo, Bunkyo-ku, Tokyo 113 8654, Japan 


Abstract: Radioactive waste disposal is important facility for human and environment in the world. Compacted bentonite in 
radioactive disposal engineer barrier design really experience hydration effort as decreasing of suction during long-time. Hydration 
effort develop macro-micro void structure in bentonite under deeply geological environment. The bentonite occurred uncertainly 
problems or translation in various experimental interaction boundary conditions such as thermal-hydration-chemical condition. To 
detect accumulation of deformation or changing of bentonite behaviour due to these processes is important that the modified 
experimental methods are required. In addition, to interpret laboratory experimental results combine to establish mathematical 
modelling in possible. The overall investigation or performance of the bentonite have contributed to represent the intrinsic properties 
of engineer barrier systems. This study focused on changing of properties of unsaturated compacted bentonite related to hydration 
effort due to increasing of relative humidity. Changing of some properties revealed to become instability or uncertainly problems in 
practice. Soil-water characteristic curve was measured with considering of various temperatures using vapor pressure technique. 
Swelling pressure and creep behaviour such as mechanical components were described with hydration effort. 


Key words: Bentonite, suction, soil-water characteristic curve, swelling pressure, creep deformation. 


1. Introduction repositories. The safety design analyses for a deep 


, , , geological repository was modified based on a 
Radioactive waste disposals have been produced , 
, en comprehensive methods or measurement program of 
from atomic plant, and decommissioning of the 
expansive materials such as bentonite, expansive clay 
nuclear power plants have supplied severe problems to 
and host rock. The general procedure has been 
environment. Deep geological disposal of radioactive n 
Ma evaluated by many laboratory tests and in-situ 
waste or some of decommissioning of the nuclear , o 
, , , investigation [1-8]. 
power plants are science, technically and societally 
The bentonite is one of interesting materials which 
challenging duty. Geotechnical or Geo-environmental 
consist of engineered barrier for deep geological 
engineers in the world have approached some kind of 
repository. There are so many experimental reports in 
mandate with the necessary honour and usefulness for 
area of physical, science or technical. Following 
highly generational protection of human. Actually, o. , 
, , , , preliminary laboratory test, the work on swelling 
mount of radioactive disposal waste with some levels 
À , ae pressure test updated the enormous database or 
arises according to the nuclear energy legislation or as 
o. , interpretations on the experiments that a number of 
result of powerful economic activity. Unsaturated soil l 
oe synthesis reports were presented [9-16]. 
mechanics is required to response within the accurate 
Working developments must advise the establish 
framework for planning, construction, operating and a 
safety program to barrier system for submission of 
maintenances of the considerable deep geological oo, l , 
license application: one low and intermediate level 


T nd one for spent fuel, high-level and long li 
Corresponding author: Tomoyoshi Nishimura, Professor, waste and one fo spent tuel, high-level a d lo 8 ved 


main research fields: unsaturated soil mechanics. intermediate-level waste. It is necessary to evaluation 
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Fig. 1 SWCC with various temperatures. 


of bentonite properties in unsaturated condition using 
suction control method or soil moisture retention 
concept [4, 17-27]. Many reports regard to soil-water 
characteristic curve or suction controlling procedure 
have been initiated to expand the mathematical 
simulation models [28-33]. Then, shear strength 
mechanics or experimental method with suction 
controlling supported the interpretation of unsaturated 
bentonite properties as scientific and technical 
background [34-46]. 

Recently, many data base obtained from geological 
field investigations in boreholes notice that there are 
uncertainly problems in the engineered barrier system 
such as hydration effort or thermal efforts. Hydration 
effort is one of factors associated to changing of 
relative humidity environment [47-53]. Thermal 
conductivity is one of parameters on 
Hydro-Thermal-Mechanical-Chemical simulations. 
hydration effort induced the 


deformation to unsaturated bentonite. The deformation 


Particularly, sever 
had closely creep behavior [54-61]. 

This study focused on hydration effort to properties 
of compacted bentonite, which was used as one 
component of radioactive disposal waste barrier 
system. Hydration induced volume deformation of 
compacted bentonite that changing of soil moisture 
cause, and vapour pressure technique was useful to 
control suction. Soil-water characteristic curve of 
compacted bentonite was indicated with both drying 


pass and wetting pass that increment of reduction in 


suction were applied to bentonite at term long-time. 
Also, swelling pressure was measured to compacted 
bentonite subjected to different suction. It was effort 
that suction decreased before swelling pressure test. 
All of these properties were measured using vapor 
pressure technique, which one of suction controlling 
methods. Also, the developed apparatus was used to 
appear creep behaviour of unsaturated compacted 
bentonite. After reduction of suction of bentonite, it 
was indicated that bentonite approached up to large 


deformation. 


2. Soil-Water Characteristic Curve 


2.1 Soil Material 


The water retention activity is important feature 
which described as the soil-water characteristic curve 
(i.e. SWCC). The SWCC was one of key parameters 
for establishing of mathematical models, and basic 
physical parameter. Mechanical and hydration 
properties were evaluated using the SWCCs in the 
general soil material. The soil-water characteristic 
curves for sodium type bentonite were measured, 
which was Kunigel V1. SWCC test was conducted out 
using vapor pressure technique. The soil-water 
characteristic curve of sodium bentonite was indicated 
in Fig. 1. The specimens had no compaction that was 
powder condition. The determined SWCCs had a wide 
range from 2.8 MPa to 296 MPa which was obtained 
using the various salt solutions. Using vapor pressure 
technique for high suction measurement is no longer 
special suction application in unsaturated soil 
mechanics. In this testing, SWCCs was controlled 
with various temperature, which had a range from 20 
degrees to 60 degrees, and difference was ten degrees. 
Obtained water contents in SWCC decreased with 
suction in drying process, and influence of 
temperatures on water content was quite small. It was 
possible to define one curve having an inflection point. 
So, unique opinion associated with increasing of 


temperature was not confirmed from these data sets. 
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The compacted unsaturated bentonite was prepared 
that changing of both water content and degree of 
saturation were measured under RH of 98%. RH 98% 
equilibrium with suction of 2.8 MPa, which was less 
than initial suction of compacted bentonite. Hydro 
effort imposed to bentonite, and soil moisture 
increased according to time. Phenomena-induced 
hydration were shown in Figs. 2 and 3. Increment of 
water content was obtained to measure changing of 
weight of specimen with time, and degree of 
saturation was calculated using changing of weight 
and changing of both diameter and height at any time. 

The compacted bentonite had a dry density of 1.600 
g/cm? before applying of hydration. Increment of both 
water content and degree of saturation proved to 
induce the deformation of bentonite by increasing soil 
moisture due to hydration effort. 

It is important to investigate the influence of salt 
solution on water retention activity of compacted 
bentonite. This study conducted out SWCC test for 
saturated bentonite in salt water. Saturated bentonites 
were prepared for measurement of SWCC that one 
was swelled in the distilled water, and another one 
was swelled in the salt water (i.e. NaCl water). Salt 
water had 3.5% in concentration, which was similar to 
nature sea water. Each bentonite was dry density of 
1.600 g/cm? at initial condition and had no mixture 
sand. Compacted bentonite with water content of 
8.0% was absorbed in each bentonite was dry density 
of 1.600 g/cm” at initial condition and had no mixture 
sand. Compacted bentonite with water content of 
8.0% was absorbed in steel mould, and initial 
specimen size was a height of 2.0 cm and a diameter 
of 6.0 cm, respectively. Absorbing a bentonite 
specimen was conducted under constant initial volume 
during over one month. After swelling, it was 
confirmed that each specimen approached to 
saturation in the steel mould. Each saturated bentonite 
was placed into glass desiccator with salt solutions. 
Variety salt solutions produced specified relative 


humidity that was vapor pressure technique. This 





testing program used seven different salt solutions, 
and controlled relative humidity had a range from 98% 
to 11%. Using relative humidity ranges corresponded 
that suction had a rage from 2.8 MPa to 286 MPa. 

As first testing step, saturated bentonite was placed 
into RH 98% environment in glass desiccator, and this 
step maintained over one month. On evaluation of 
equilibrium with relative humidity to bentonite 
specimen, a period of one month was in generally 
accepted in SWCC test procedures, was recognized as 
one of practice methods in conventional. When 
specimen took an equilibrium, some parameters were 
measured which were weight, diameter and height for 
specimen, and determine water content and degree of 
saturation to controlled suction. Step by step, suction 
was decreased up to 286 MPa due to replace salt 
solutions in drying process. Relative humidity in glass 
desiccator according to be reduction that relative 


humidity was 11% at end of drying process. All of 
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Fig. 3 Increment of degree of saturation under according 
to hydration. 
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suctions were seven different values in drying process. 

Subsequently, wetting process-induced by reduction 
of suction was commencement to all bentonite 
specimens experienced drying process. Suctions on 
process in wetting procedure were quite same that it 
was easy to comparison to observing test results at 
drying process. Above mentioned, physical properties 
such as water content and degree of saturation were 
determined with same testing methods. 

Both water content and degree of saturation such as 
parameters were used to describe each SWCC of 
bentonite (Fig .2 from 4 to 7). On drying process, soil 
moisture was obviously large for bentonite including 
salt component up to suction of 9.8 MPa. Beyond 
suction was 9.8 MPa, reduction of water content with 
increment of suction was indicated, and thus results 
were similar between bentonites with and without salt 


water at drying process. Also, water contents were 
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Fig.5 Changing of water content (salt water). 
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Fig. 7 Changing of degree of saturation (salt water). 


measured near zero value regardless of imposing salt 
components at suction of 296 MPa. 

On wetting process, all bentonite specimens 
described to increase soil moisture regardless of 
experience with and without salt water that it was 
widely accepted in water retention mechanics in 
unsaturated soils. Water content and degrees of 
saturation can seem to be not same with those data 
sets of bentonites in drying process that exist much 
hysteresis in hydro-processing. Particularly, the 


bentonite with salt component indicated large 
increment of soil moisture with reduction of suction 
comparison with observing results in de-stilled water. 
As results, few hysteresis was observed in the 
bentonite having salt components that strong efforts 
contributed to water retention activity associated with 


bentonite properties. 
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3. Swelling Pressure of Unsaturated 


Bentonite Subjected to RH 


Swelling pressure is one of important key 
parameters for construction of radioactive waste 
disposal system which is supplied saturation process 
under deep ground during long term period. Many 
reports regard to swelling properties of variety 
bentonites had been established that almost of 
swelling tests were conducted to supply water 
absorbing to unsaturated bentonite directly. It is sure 
that barrier structure constructed using unsaturated 
bentonite was initially unsaturated condition, and 
considerable slow water flow is predicted as results of 
many mathematical simulation models. To reach to 
completely saturation situation spend extremely long 
times. While unsaturated bentonite become toward to 
saturation, hydration effort make slowly expansion in 
changing of suction of bentonite under deep tunnel. It 
was no consider that deformation of unsaturated 
bentonite occurs according to suction reduction (i.e. 
increment of relative humidity). 

A swelling pressure testing apparatus was 
developed to determine swelling pressure in a constant 
relative humidity environment. The apparatus consists 


of a triaxial chamber and relative humidity control 


circulation system. The modified swelling pressure 
testing apparatus was shown in Fig. 8 that consisted 
mainly of a triaxial chamber, a pedestal, a steel mold, 
a double glass burette, a differential pressure 
transducer, a difference displacement sensor, load cell 
sensor and relative humidity control circulation 
system. The relative humidity control circulation 
system was established using a conventional pump, 
along with a small chamber with salt solution. The air 
flow maintained a constant relative humidity 
surrounding the compacted bentonite. All compacted 
bentonite specimens were placed into steel mould. 
Moving water in the double glass burette due to 
absorption in bentonite was permitted at the low 
portion of specimen. Absorption water volume change 
in the double burette was measured using the 
differential pressure transducer. Initial total volume 
was maintained constant for all compacted bentonite 
specimens. 

Sodium bentonite was used for measurement of 
swelling pressure. Silica sand had uniformity grain 
size distribution which was mixed with bentonite at a 
ratio of 30% by dry weight. The specimen was 
statically compacted in stiffness steel mold at an initial 


water content of 5.9%. The compacted bentonite 
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Fig.8 Modified swelling pressure test apparatus. 
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specimen had a dry density of 1.600 g/cm’ as target 
value. The height of specimen was 25.5 mm. 
Compacted bentonite specimens with two different 
soil suctions were prepared for swelling pressure tests. 

Initial compacted bentonite had a soil suction of 
105 MPa. Soil suction of 2.8 MPa was imposed using 
vapor pressure technique and used salt solution was 
Potassium Sulfate. 

Soil suction of 2.8 MPa corresponds to relative 
humidity of 98%. The gravimetric water content 
increased that vapour moisture in RH 98 was absorbed 
into void structure of bentonite having suction of 105 
MPa. As result, the gravimetric water content of the 
compacted bentonite sample subjected to relative 
humidity of 98% increased to 9.7% from 5.9%. Also, 
deformation in expansion occurred, grew with time. 
Finally, deformation of bentonite specimen took large 
expansion in vertical direction, and was 9.5% in 
swelling. During hydration process, void structure 
constructed macro-micro structure produce complexly 
expansion in clay structures together. Changing of 
total volume was useful to determine dry density. The 
dry density of the sample after equilibrium soil 
suction of 2.8 MPa decreased to 1.427 g/cm? from 
1.600 g/cm’. 

Swelling pressure of bentonite with two different 
suctions were measured under constant volume, which 
changing of a height of specimen was not permitted 
through swelling test. 

Fig. 9 shows observing the swelling pressures with 
water absorb for compacted bentonite at initial water 
content of 5.9%, and specimen had an initial soil 
suction of 105 MPa at beginning test. The measured 
swelling pressure described rapidly growing during 
almost 50 hours periods. Strong increment of swelling 
pressure seems to be one of characters on swelling 
properties for compacted bentonite with high suction. 
Determined swelling pressure at 50 hours was 220 
kPa. Behind 220 kPa in swelling pressure, bentonite 
described traditionally small reduction with time. 


Reduction continued smoothly till elapsed time was 


350 h, and measured swelling pressures had date sets 
with both increment and decrement such as a 
frequency cyclic motion. Once again, the bentonite 
showed small increment, which appeared increment 
and decrement in swelling pressure. Its smoothly 
increasing was maintained till swelled bentonite took 
equilibrium, and elapsed time was nearly 1,000 hours. 
Confirmed maximum swelling pressure corrected as 
267 kPa. 

The compacted bentonite with soil suction of 2.8 
MPa shows smoothly increasing of swelling pressure 
at beginning of test as shown in Fig. 10. There was 
quite difference comparison to swelling behaviour of 


initial suction of 105 MPa. It was obviously that 
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Fig. 10 Processing of swelling pressure (hydration). 


98 


performed swelling pressure was smaller than that of 
suction 105 MPa at beginning that was revealed as the 
influence of suction of bentonite before absorbing. 
Growing of swelling pressure was extremely slow till 
elapsed time was 350 h that was entirely difference 
against to swelling pressure increment mentioned in 
Fig. 9. 

After the elapsed time of 350 hours, accumulation 
of swelling performance became to be developed, and 
it reached to 200 kPa at 650 hours. Beyond 750 hours, 
bentonite seemed to be equilibrium that 206 kPa was 
measured as maximum swelling pressure. Bentonite 
showed that hydration effort caused reduction of 
swelling pressure, and at same time growing processes 
of swelling pressure was quite difference between 
with and without hydration effort before swelling 
process. Consequently, degree of saturation obtained 
from water content at end of test reached to saturation 
regardless of magnitude of suction value before 


swelling. 


4. Creep Test 


This study mad above mention that hydration effort 
such as reduction of suction improved both water 
retention capacity and swelling performance in 
compacted bentonite. Revised basic properties to 


compacted bentonite was induced at macro-micro 
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structure as extremely fine structures. This testing 
conducted out creep test and finding relationship 
between main chief factor in hydration effort and 
strength resistance in mechanical properties. It is 
defined in generally for the soil mechanics that creep 
test is conducted under constant effective stress to 
saturated soils. The changing of macro-micro void 
structure due to hydration effort had closely the 
deformation on creep behavior. The modified creep 
test apparatus was shown in Fig. 11. The apparatus 
employed a conventional cyclic relative humidity 
control system, which was possible to apply a required 
RH using vapor pressure technique. This apparatus 
consists of triaxial chamber, air cyclic flow pump. A 
dynamic activity of the conventional pump had at least 
10 kPa pressure, which maintained steady air flow 
through the system. 

Creep stresses were determined using unconfined 
compressive strength in experimental data sets what was 
deviator stress without confining pressure. Controlled 
relative humidity was either 98%. In case of RH 98%, 
suction of 2.8 MPa corresponding to RH 98% was 
considerable lows level suction comparison with 
suction of bentonite in initial condition. Used specimen 
had both a height of 100 mm and a diameter of 50 mm, 
respectively. Also, dry density was 1.600 g/cm? as 


target value. 
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Fig. 11 Modified creep test apparatus. 
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Creep stress 
(deviator stress): 
182.8 kPa 


Compaction condition 
Water content: 8.0 %0 
Dry density: 1.600 g/cm? 


Shinkage 





Fig. 12 Case of no hydration for creep test. 
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Fig. 13 Case of application of hydration for creep test. 


Observing axial strain under subjection of external 
vertical loading was shown in Fig. 12. Prepared creep 
stress was 182.8 kPa and was maintained till end of 
test. Then, positive value in axial strain expressed that 
specimen occurred compression deformation. The 
axial strain in shrinkage was approximately 1.0% at 
beginning of external loading, which was produced by 
deviator stress as mechanical action. The axial strain 
was remained a period of 30 days (i.e. up to end of 
test). It was able to predicate that the specimen 
described growing shrinkage deformation under no 
hydration environment. 

As other case, axial strains with elapsed time was 
indicated as shown in Fig. 13, and applied creep stress 
was 364.5 kPa more than that mentioned in Fig. 12. 
Though bentonite had creep stress of 364.5 kPa under 
no hydration till fifty days, compression deformations 


were measured at commencement of supplying 


vertical load. While fifty days, axial strain described a 
stable, and the strains were similar with measured 
strain on beginning of loading of stress. There were 
not observing increment of shrinkage strain such as 
creep deformation. 

Subsequently, the hydration was applied to 
bentonite due to air flow having RH 98%, and suction 
surrounding of specimen approached rapidly 2.8 MPa 
in suction. As point of view regard to presence or 
absence in hydration effort, beyond fifty days the 
bentonite occurred the expansion distinctly with days. 
It was sure that the effort of hydration was revealed 
with measured changing of vertical strain. The 
processing of hydration effort developed macro-micro 
structures in bentonite, and large increment of 
expansion was described from forty days to fifty days. 
At fifty days, 


comparison with 


specimen had large expansion 
height that 


approximately 


measured 
1.0%. 


Subsequently, the bentonite evidenced the shrinkage 


initial 
expansion strain was 
deformation such as up to destruction. Thus, axial 
strain data sets indicated straight line in vertical 
direction. Crushed specimen was clearly observed in 
chamber. Then, the hydration effort can be defined as 
decrement of suction in stress variables that induced 
the destruction combined large deformations for 


unsaturated bentonite. 
5. Conclusions 


This study represented properties of bentonite such 
as extremely expansive soils on experimental resulting. 
The hydration effort to compacted bentonite was 
referred as one of significant impact factors, which 
induced uncertainly conditions on engineered 
properties to unsaturated bentonite. On supporting 
mechanics, experimental 


unsaturated soil some 


laboratory tests were conducted used modified 
apparatus that were SWCC test, swelling pressure test 
and creep test for unsaturated compacted bentonite. 
These tests had enough advantage to investigation 


above mentioned physical and mechanical properties 
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of compacted bentonite. Particularly, interpretation 
incorporated suction controlling methods were useful 
to understand the influence of hydration effort such as 
reduction suction on bentonite properties: 

(1) The bentonite had high water retention activity 
at high suction ranges such as larger than 2.8 MPa in 
suction. The hydration effort induced expansion of the 
volume under reduction of suction, and, increment of 
soil moisture occurred at same time. Also, reduction 
of suction was controlled using vapor pressure 
that 
controlled in mythology. 


technique relative humidity was directly 

(2) The obtained soil-water characteristic curves of 
bentonite can affirm that the influence of temperature 
on water retention activity was negligible. Other hands, 
the salt water produced to be high soil retention 
activity that bentonite absorbed, swelled in the salt 
water. 

(3) Due to suction reduction, the swelling pressure 
with time indicated smoothly growing comparison to 
initial condition (i.e. before applying of hydration). As 
results, peak swelling pressure decreased that the 
transmutation of macro-micro void structure in 
bentonite immediately related to performance of 
swelling behavior. 

(4) Through creep test, when hydration effort such 
as suction reduction was applied to compacted 
bentonite, it was obviously that axial strain in 
expansion direction was observed, and unsaturated 
compacted bentonite reached up to failure as result of 
maintain of both hydration effort and creep stress 


loading. 
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Abstract: In this paper, an overview of an important feature in statistics field has shown: the stepwise multiple linear regression. 


Likewise, a link between stepwise multiple linear regression and earthquakes localization has been descripted. Precisely, the aim of 


this research is showing how stepwise multiple linear regression contributes to solution of earthquakes localization, describing its 


conditions of use in HYPO71PC, a software devoted to computation of seismic sources’ collocation. This aim is reached treating a 


concrete case, that is computation of earthquakes localization happening on Mount Vesuvius, Italy. 
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1. Introduction 


Linear regression is probably the simplest and, at 
same time, the most well-known method of research 
of a mathematical relation between an explanatory 
variable and a dependent one of physical quantities. A 
particular specification is about the case in which 
explanatory variables are more than ones, that is the 
case of multiple linear regression. Basic concept is to 
establish a priori that physical quantities which are at 


stake are linked with each other in a linear way by 


means of a mathematical relation, named linear model. 


The goal is finding the parameters of this model, 
starting from observed data. A particular type of these 
statistical models is the so-called SMLR (Stepwise 
Multiple Linear Regression). In this regression 
typology, computation of model’s parameters happens 
by means of an algorithm. In this paper, a description 
of F-test as a good criterion for implementation of this 
linear regression typology has been shown. Otherwise, 
use of SMLR to solve an important physical problem, 
earthquakes localization, has been shown. Indeed, to 
establish in a exact way geographical coordinates of a 


seismic event is very relevant to obtain all the 
Corresponding author: Giuseppe Pucciarelli, Ph.D., 
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essential information about it. In particular, an 
essential point is finding geographical coordinates of 
seismic events which happen on volcanic areas. If 
these events happen in a tight period of time and their 
hypocenters are collocated along volcanic pipe, it 
could possibly be an imminent magma eruption. This 
is the reason because an application of computing 
earthquakes localization in Mount Vesuvius by means 
of SMLR has been shown. Mount Vesuvius is an 
active volcano which is placed about 9 km east of 
Naples, Italy: since a population of 3,000,000 people 
live in its foothills, it 1s one of more monitoring 
volcanoes in the world. The entire process to locate 
specific earthquakes which occurred along Mount 
Vesuvius crater by means of SMLR technique (or 
better yet, by means of a software named HYPO71PC 
based on SMLR technique) has been described. In this 
paper, this kind of explanation is preceded by a 
generic overview about earthquakes localization 
problem and it is followed by a concise description of 
the applications which have earthquakes localization 


as starting point. 


2. Statistical Model: The Stepwise Regression 


Mathematically speaking, the regression is a 


statistical process which allows obtaining an estimate 
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of an assumed value from a dependent variable Y as a 
sist Xn. The 


most common type of regression is the linear one. 


function of n independent variables X; , 


This is based on the idea that for each sample of 
observations, there has been a determination of the 
variable Y and n determinations of the variables X; , ..., 
Xn. The goal is to find a linear relation between these 
variables. As is well-known, the most simple linear 
relation is that given by an affine transformation linear 
that has the Eq. (1): 

Y= Bot B X,+ E; (1) 
where i varies between the observations, 1.e. i = 1, ...., 
n; Y; are the dependent variables; X; are the 
independent variables, called regressors; fy + 6X; 1s 
the regression line, called also regression function; €; 
represents the statistical error. 

Eq. (1) has the Eq. (2): 

Y=B,+B,X+e (2) 
in which the vectors Y, X and € are n-dimensional. The 
values fo and pı represent the parameters of the model 
which have to be estimated by the obtained 
observation samples (so, they are deterministic); 
instead, the vector € contains only stochastic 
parameters to be determined, because for every 
assumed value from X, it exists in an entire probability 
distribution that contains values of Y (for this reason, 
it is impossible to know with certainty the value of Y). 
Thus, the relation Eq. (2) gives a random variable, 
where its probability distribution and characteristics 
are determined from the values of X and the 
probability distribution of €. Among all of the methods 
of linear regression known, the most used is surely the 
method of least squares. 

The type of linear regression as soon as view is said 
simple linear regression, because it is available to 
have only one independent variable. If instead, there 
are more independent variables that can be used for 
the determination of the dependent variable Y, we talk 
about multiple linear regression. In this case, the 


equation for the previous problem, becomes: 
Y,=B t+ B, Xut B Xat t BX yt 8) 


where with respect to Eq. (1), the only change is that 
the regression line is Po + PiXii + p2Xhi + ... + PkXki + &i. 
This equation has a matrix form: 

Y=XBteE (4) 
in which X is a k x n matrix, while the vectors Y, € are 
n-dimensional and f is k-dimensional. A type of 
multiple linear regression very used for the 
applications in the last years is the so-called stepwise 
regression [1], discovered from Draper and Smith in 
1966 [2]. It is especially used when, in the process of 
regression, problem is solved by means of 
independent variables which fail to give information 
on the variance of the dependent variable: stepwise 
regression is able to solve this sort of problems. In this 
method, the total variance is partitioned into two parts: 
the first one which can be explained by the 
independent variables (explained variance) and the 
second one which cannot be explained by them 
(unexplained variance) [3, 4]. In order to do this, 
5 An, deemed 


suitable to the explanation of the variance of the 


given n indipendent variables X, 


dependent variable, it is possible the extraction step by 
step (hence the name) of the independent variable that 
has the highest percentage of the explained variance. 
Now considering the step k + 1: the variabile X, +1 1s 
related with X, of the previous step: thus, Yx is called 
residual variable, because it is given by Y minus the 
linear combination of the independent variables 
already present in the model. In practice: 
PaSa Bo- BX i7. BX; (5) 
The variable X;,.; introduced in the model will be 
the one that will show the highest correlation 
coefficient, in absolute value, with the residual 
variable. After establishing the criterion for choosing, 
an algorithm chooses the independent variables 
suitable for the explanation of the variance of the 
dependent variable, after an initial setting of a 
threshold value with which comparing the diffierent 
coefficients [1l, 2]. 


procedure is not completely reliable (even considering 


correlation However, this 


that the threshold value is subjective), therefore the 
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diffierent correlation coefficients have to satisfy the 
following null hypothesis: 

Hy: 6,+ 1 (6) 
where b,+1 is the so-called regression coefficient, 1.e. 
the slope of a line obtained using linear least squares 
fitting [5, 6]. So, Eq. (6) affirms that the regression 
coefficient of the independent variable Xk+ to input in 
the model is zero, after the introductions of Xj,..., Xx. 
This hypothesis is verified through the F value of 
Fisher, given by 


R? 

MSR_ k 
f= a 
m-— k- 1 


in which MSR and MSE are respectively the mean 
square due to regression and residual; furthermore, the 
coefficient of determination R? is given from: 

j VAR 


R = ]- res 
VAR, (8) 





where the numerator (respectively denominator) is the 


sample variance of the estimated residuals 
(respectively dependent variable). It is clear that 0 < 


R° <1 [7]. 


3. Limits of the Model and Implementation 
of the Algorithm 


The F value of Fischer is based on a single 
hypothesis: in the stepwise regression, it is hugely 
violated [7]. As a result, models are too complex and 
other parameters could be biased towards 0. Therefore, 
diffierent alternatives for the stepwise regression are 
developed. The following methods perform better than 
stepwise, but their use is not appropriate for statistical 
knowledge [8]: the most important are LASSO, LAR, 
CROSS VALIDATION and 


procedures (not yet enough reliable). Even though no 


some experimental 
one methods could substitute for statistical experience, 
LASSO and LAR (especially) are much better than 
stepwise: these two methods allow us to consider a 


wide range of options [8]. Through the F value of 


Fischer, the model with the introduction of the 


independent variable X4; is still statistically 
significant compared to the previous one, which 
.., Ak: this 


value is also called F to enter. The introduction of new 


included the independent variables X4, 


variables could, however, result in an information loss 
for some variables introduced previously: to check 
this phenomenon, authors compare the correlation 
coefficient with another threshold value, known as F 
to remove: if they are present, these are immediately 
deleted. Usually, almost all programs have the 
possibility to follow step by step the regression 
procedure, through information about the variables 
and the signalling of F to enter and F to remove which 
the program utilizes for this particular step. A further 
option frequently present in the various programs 
concerns the possibility of choice, on the part of the 
user, of variables to be introduced or delete in the 
model: this choice is made through the control of the 


values of F to enter and F to remove. 
4. Earthquakes: The Geiger’s Method 


Earthquakes’ localization constitutes a particular 
scientific problem where the goal is computing the 
coordinates of a point situated inside the Earth from 
which seismic waves are generated. This point is 
known as the hypocenter of an earthquake and its 
outside projection is named epicenter. Hypocentral 
coordinates which authors want to ascertain are four: 
three are spatial (latitude, longitude and depth) and the 
other one is the origin time of seismic waves 
generation. Historically, earthquakes localization was 
carried out by means of analysis of outside 
consequences produced by a seismic event. Now, 
thanks to seismic network development and a more 
detailed knowledge of interior structure of the Earth, 
computing of the four earthquakes localization 
parameters could be done in a more meticulous way. 
Earthquakes’ localization is based on classical linear 
motion law, that is: 


x= vt (9) 
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where: 

e x is the hypocentral distance, that is the distance 
between hypocenter and seismic station; 

e y represents the seismic waves velocity; 

e tis the seismic waves traveltime. 

Describing the process in a more practical way, 
earthquakes localization needs three important 
features, 1.e.: 

e knowledge of seismic stations’ geographical 
coordinates; 

e picking of seismic phases on seismograms and 
computing of their traveltimes; 

e knowledge of velocity structure between 
hypocenter and seismic stations. 

Therefore, earthquakes’ localization is a typical 
inverse problem, that is a problem in which solution 
consists in an evaluation of parameters of a 
determinate model (they must be in good agreement 
with observations). In this particular situation, 
parameters are the hypocentral coordinates, while the 
observations are picked seismic phases and model is 
an Earth hint between hypocenter and seismic stations. 
Rewriting Eq. (9) for seismic phases P and S, 
respectively, that is the first waves generated by a 
hypocenter which arrive to a seismic station, 
following equations are obtained: 

A= v (t ty); A= v(t; to) (10) 
where A is the hypocentral distance, vp and vs are the 
velocities of P wave and S wave respectively, tp and ts 
are the P traveltime and the S traveltime and fo is the 
seismic event origin time. Manipulating in an 


opportune way Eq. (10), Eq. (11) is obtained: 
KA (11) 


Replacing Eq. (10) in Eq. (11), Eq. (12) has been 


obtained: 
Vp Vp 
t— t, = |—- 1t, ,—-|{—-1liz 
S FP LS. tp E to (12) 


Eq. (12) can be understood as an equation of a line 


in a cartesian plane where fpis the abscissa and (ts — tp) 
is the ordinate. In this way, through a simple linear 
interpolation, parameters can be found. By means of 
slope, it is possible to compute vp/vs , while the 
y-intercept allows the computation of seismic event 
origin time tọ. From values of these specific 
parameters, it is easy to obtain A through Eq. (11). 
Earthquakes’ localization can be realised in a graphic 
way by means of so-called circle method [9]. If at 
least three seismic stations are at our disposal, for each 
station the epicentral distance has been computed. 
Then, a circle with the centre arranged on a seismic 
station and with radius equivalent to epicentral 
distance is drawn. The point of intersection among all 
the circles is the localization of an earthquake. A 
modification of this method is the following. If only P 
seismic phases and three stations are computed, circle 
method is applied in the same way as previously 
described for two of three stations. In this case, 
earthquakes’ localization is the centre of a circle 
tangent at two drawn circles and passing through the 
station. Both 


interpolation and circle method show us the following 


third remaining seismic linear 
set of problems: 

e Earth between hypocenter and seismic stations is 
considered as a homogeneus halfspace where the 
velocity of seismic waves is constant. Obviously, this 
is not the real case. 

e Earth’s depth is completely disregarded. 

e Consideration of a spherical Earth is completely 
disregarded. 

In 1912, L. Geiger introduced an algorithm to solve 
the problem of earthquake localization, considering 
three criticisms mentioned before. This method, justly 
called Geiger’s method, is still used as basis of 
software for earthquakes localization [10]. Starting 
point of this method is model of representing Earth as 
a sequence of parallel layers where in each one of 
them, seismic waves velocity is constant. Starting 
from a trial solution for localization, it is named mo(Xxo, 


Yo, Zo, to). After, there is the procedure of computation 
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of the so-called residuals. For the i-th seismic station, 
Eq. (13) is defined: 

R= do, (13) 
where R; is the i-th residual that is the simple 
diffierence between theoretical traveltime T; and 
observed (by means of seismograms interpretation) 
traveltime ti. Writing T; in function of mo, Eq. (13) 


becomes: 
R=T + l E A EE oe eee o E 
p zio Xi) + (Vom a a a 
(14) 
where v is the velocity of considered seismic wave. 


Now, rewrite Eq. (14) and apply on it an expansion 


in Taylor series at mo. In this way, Eq. (14) becomes: 
R,~R,(m,)+ a,Axt+b,Ayt+c,Azt+At (15) 


OR, aR OR. 


where a; = a b; = Dann Ci = Jz, it all 
partial derivatives calculated in mp and Ax = x — xo, Ay 
= y — Vo, AZ = Z — Zo, At = t — tọ. 

Finally, a least squares method to make residual 
minimum is used. Therefore, hypocentral parameters 

n 
are computed when » R? is minimum with n 
i=] 
seismic stations. 

Summarising, to localize an earthquake, this 
procedure should be: 

(1) Acquisition of seismic phases and of observed 
traveltimes by means of seismograms. 

(2) Geographical coordinates of seismic station 
which have minimum observed traveltime become 
spatial coordinates of trial solution mọ. While authors 
assume as fo, the P seismic wave traveltime registered 
at that seismic station deprived of a quantity 
depending on velocity model authors choose. 

(3) Computation of residuals. 

(4) Appropriate correction of trial solution. 


(5) New computation of residuals. 
(6) Repetition of (4) and (5) until when » R° is 
=I 


minimum. 
It is appropriate to emphasize how solution’s 


goodness depends on an adequate seismic rays 
coverage, a sufficient knowledge of velocity model 


and a careful picking of seismic phases. 


5. Application: Examples of Earthquakes 
Localization on Mount Vesuvius 


Mount Vesuvius is an active stratovolcano which is 
situated about 9 km east of Naples, Italy, and a short 
distance from the shore. Since a population of 
3,000,000 people live in its vicinity, it is one of more 
monitoring volcanoes in the world. And one of 
monitoring strategies is earthquakes localization. 
Indeed, earthquakes localization is very important in 
volcanic areas because a contingent high 
concentration of seismic events could be due to 
magma rising along the volcanic pipe. Magma rising 
causes a gradual deformation of volcanic pipe rocks 
which can trigger seismic events. And often magma 
rising could forerun an eruption. Therefore, accuracy 
of earthquakes localization in volcanic areas could 
prove to be fundamental. In this paper, a computed 
earthquakes localization of 2,315 seismic events 
which occurred on Mount Vesuvius in a period from 
06/12/1986 to 31/12/1999 and was registered by 41 
seismic stations belonging to OV-INGV Mount 
Vesuvius monitoring network is shown. It is important 
to underline that only P and S seismic phases have 
been chosen for localization. Starting point is choice 
of an opportune velocity model to schematise Earth. 
This last one is a standard 1-D velocity model for 
Mount Vesuvius area (see Fig. 1) which has already 
been in other previous works with satisfactory results. 
This model schematises Earth in a depth interval 
between 0 and 30 kilometers. Model is composed of 
six layers and in each of them P wave velocity is 
Then, 
localization 1s computed by means of HYPO71PC 
software. This software is a FORTRAN code and it 
represents one of more sofisticated versions of original 
code HYPO71PC. Other similar 
HYPOELLIPSE and HYPOINVERSE. This software 


assumed to be constant. earthquakes 


versions are 
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is nowadays widely used for earthquakes localization. 
The reason for this is the following one. This software 
makes use of stepwise multiple linear regression to 
carry out six points of earthquakes localization 
procedure which have been previously described in 
Chapter 5. In particular, it carries out the point (6) and 
it shows lots of benefits in accuracy and convergence 
of solution. In this particular case, an input file named 
VESUVIO.PHA is supplied to Hypo71PC. This file is 
organised in the following way: First of all, it presents 
the command HEAD followed by the number of 
characters used to write a heading above each 
earthquake in the output. This command is put equal 
to 20. Then, the command RESET 1s written. This 
command allows setting values of the software test 
variables in a way which is straightly connected to 
input data. In this particular case, seven test variables 
are test. Among these, the critical F value for stepwise 
multiple regression is set. In Hypo71PC, this value 
should be included between 0.5 and 2 according to 
number and quality of available P and S traveltimes. 


Indeed, if this value is lesser than 0.5, then 


coefficients matrix of Eq. (15) could be ill-conditioned. 


Instead, if this value is greater than 2, then Geiger’s 
iteration may finish prematurely and therefore a good 
localization could be not obtained. For avoiding it, 
critical F value is set equal to 0.5 because there are not 
lots of available S traveltimes. Then, geographical 
coordinates (latitude and longitude) of seismic stations 
which have recorded earthquakes to localize must be 
written. After the realization of this procedure, 
velocity model must be choosen. Hypo71PC procures 
two options: station delay model and variable first 
layer model. In first case, there is the simple addition 
or subtraction to traveltimes of a number which is 
stated for each seismic station. In second one, this 
stated quantity—previously described—is interpreted 
as a thickness to add to first layer. Therefore, each 
seismic station has the same P-wave velocity but 
thickness of first two layers under each seismic station 


varies. This typology of model is often chosen when 


there is the case of different travel paths, that is to say, 
paths do not come from an uniform reason. 
Fortunately, this is not the case and therefore station 
delay model is chosen. Following step is writing a 
record in which parameters of trial solution and 
control indicators for the software are indicated. In 
these parameters of the trial focal depth, the distance 
from epicenter and others concerning the earthquakes 
localizations are set (see Ref. [11]). It is important to 
underline how there are other possible options that can 
be activated, especially, options concerning 
geographical coordinates of a trial solution, but this is 
not the case. Therefore, for each earthquake, 
Hypo71PC chooses as geographical coordinates of 
trial solution those of seismic station which is nearest 
to earthquake localization. Finally, in the last part of 
input file traveltimes are found. They are organised in 
the following way. For each earthquake, this is the 
structure of the script: the acronym of seismic station 
which has recorded that earthquake, the symbol of 
analysed seismic phase (that is, P), P-traveltime at that 
station with accuracy down to centiseconds and, when 
they have been recorded, S-traveltime at that station 
with accuracy down to centiseconds. When all 
traveltimes records regarding to a specific earthquake 
have been written, then the character “10” in the 
record immediately following is written. This 
particular record is considered as a pause between 
traveltimes of a seismic event and traveltimes of 
another one. The output file is named VESUVIO.PRT: 
it is divided into two parts, where in the first one the 
most important thing is given from the root mean 


square error of time residuals, that 1s: 


= Vo (16) 


where R; is the quantity described in Eq. (14), while 
NO is the number of station readings used in locating 
the earthquake. The second part of output file is about 
the characteristics of seismic station: for further 


information (see Ref. [11]). The result of earthquakes 
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Fig. 1 1D velocity model for Mount Vesuvius area used as input for computing earthquakes localizations. 
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Fig.2 Earthquake localization in Mount Vesuvius area: the hypocenters are represented by these white circles in the figure. 


localization in Mount Vesuvius area is shown in 
Fig. 2. 


6. Conclusions 


In this paper, an interesting characteristic of an 
important topic in statistical field is described: the 
SMLR. Especially, the aim is to underline a particular 
link between statistics and physics, explaining how 


SMLR can be used in the solution of a very important 


physical problem: the earthquakes localization. Indeed, 
SMLR is used to implement Geiger’s method [10] by 
means of a subroutine of software HYPO71PC [11], 
12]. Its mode operation has been described through a 
concrete example as computation of some Mount 
Vesuvius earthquakes localization. Understanding how 
SMLR could be connected with this relevant physical 
problem is very important because a good earthquakes 


localization represents a starting point for other 
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geophysical research techniques as seismic 
[13, 14], 


implementation [15], good resolution of epicentral 


tomography focal mechanism 


distance [16] and obtaining of graphical images by 
means of GMT (Generic Mapping Tools) [17]. Of 
course, other methods for locating earthquakes exist, 
for example, Metropolis-Gibbs Method [18]. It is 
based on a 3D velocity model. This kind of method 
has already been used to locate earthquakes on Mount 
Vesuvius [19]. Method which has been described in 
this paper where SMLR is implemented is based on a 
1D velocity model. But 3D velocity model tests have 
already been carried out [20]. For future studies, 
development of these 3D velocity model tests using 
data treated in this paper could be appropriate. 

This map is given by the geographical coordinates, 


1.e. latitude and longitude. 
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Abstract: In Kuwait, wastewater management has gained extra attention in recent years and becomes crucial for sustainable 
industrial development sector. Among the food industry sector, dairy processing plants generate huge amount of wastewater, which is 
heavily loaded with organic and other toxic compounds. Disposal of dairy wastewater effluent without sufficient treatment can 
contaminate aquatic ecosystems. Cost efficient treatment processes that are effective in removing organic load and other 
contaminants are essential to meet stringent environmental regulations applied in Kuwait. A research study was carried out at the 
KISR (Kuwait Institute for Scientific Research) to assess the technical viability and economic feasibility of combined microfiltration 
and biological treatment system. This work presents the economic evaluation of the adopted treatment system. The results show that 
the cost of the integrated system for large scale is estimated to be US$ 1.575/m°, which is 25% less than the cost of wastewater 
transportation and treatment in conventional sewage plants. 


Key words: Wastewater, treatment, dairy, cost. 


1. Introduction environmental legislations is becoming a high priority 


a for the state of Kuwait. Dairy industrial sectors have to 
Dairy industry is very important to state of Kuwait 
n comply with the KEPA (Kuwait Environment Public 
and its products are wildly consumed by Kuwaiti , l l 
E , , , Authority) regulations for wastewater discharge and 
people. Similar to most of food industries, the dairy , , 
i Ean , reuse [5]. To satisfy these regulations the effluent must 
industry characteristically requires very large , 
E i be collected and treated to meet the quality standards 
quantities of freshwater and generates large quantites , 
set by KEPA. A research project was conducted to 
of wastewater [1]. Most of the wastewater volume , a , DSSS 
, o , assess the technical viability and economic feasibility 
generated in the dairy industry results from cleaning of , , , 
, , , of implementing MBR (Membrane Bioreactor System) 
transport lines and equipment between production , 
, , i , to treat dairy wastewater effluent generated by one of 
cycles, cleaning of tank trucks, washing of milk silos , o , , 
, , , the largest dairy companies in Kuwait [6]. This work 
and equipment malfunctions or operational errors [2, 3]. , , , 
, , , describes the experimental set-up and economic 
Disposing untreated dairy wastewater can lead to , , , 
, l l assessment of treating dairy processing wastewater by 
adverse public health and environmental impacts. , 
, , integrated MBR system. 
Examples of environmental concerns include 


objectionable odors and fly infestations that have 2. Process Description 


resulted from the disposal of the untreated effluent in , , 
The complete system for dairy processing wastewater 
open land [4]. Recently, the enforcement of , , , 
treatment involves the integration of membrane 
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treatment, environmental engineering and seawater , , 
desalination. the submerged membrane microfiltration system, 


processes. The membrane separation process includes 
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which is a polypropylene membrane filter to remove 
particles greater than 0.2 um in size from a feed 
stream. The main treatment steps in the integrated 
system are shown in Fig. 1, in which feed water 
passes through an aeration tank for biological 
treatment. An air compressor injects adequate air into 
the aeration tank. The aerobic system includes the 
process air blowers, which are installed adjacent to the 
system. The required process airflow is 227 m°/h 
introduced at the bottom of the aerobic tank through 
air scour distribution header pipes. After passing 
through the upstream flow, the mixed liquor is 
transferred by overflow to a suitable buffer flow tank 
and then pressurized to an operating pressure in 


accordance with the membrane’s design. 


3. Cost Analysis 


In KD-Cow (Kuwait Dairy Company) factory the 
dairy production process requires huge quantities of 
freshwater which is used mainly for cleaning purposes 
and, as a result, large amounts of wastewater are 
generated. Due to lack of treatment system in the 
factory, this wastewater is transported and distributed 
to remote sites through sewage tankers. The company 
uses 17 sewage tankers per day with a capacity of 30 
m°>. The purpose of this study is to evaluate the 
economic feasibility for the treatment of dairy 
processing wastewater effluent for reuse. Specifically, 
the objectives of this task are to: 

e To estimate the operating cost and total cost of 
the dairy wastewater treatment; 

e To calculate the unit costs of the treatment of 
dairy wastewater; 

e Use economies of scale to figure out the 


feasibility of the project. 
3.1 Cost Estimation 


Table 1 shows experimental plant operating data for 
the integrated system as obtained during the project. 
Table 2 provides the estimated cost details for 


operation expenses. 


3.1.1 Capital Costs 

Costs incurred on the purchase of land, buildings, 
construction and equipment to be used in the 
production of goods or the rendering of services are 
categorized under capital cost. In other words, it is the 
total cost needed to bring a project to a commercially 
operable status. Capital costs do not include labor 
costs except for the labor used for construction. 
Unlike operating costs, capital costs are one-time 
expenses, although payment may be spread out over 
many years. Capital costs are fixed and are therefore 
independent of the level of output. 

3.1.2 Operation Costs 

Total cost needed in daily operations is categorized 
under operation costs. Operation costs are variable 
costs, and are dependent on the level of output. 

Table 3 shows the total cost calculations for the dairy 
processing wastewater treatment and expenses details 
for the operation costs. The total cost is US$ 171,975. 


4. Results and Discussion 


The estimates of unit costs of the experimental 
plant for wastewater treatment are provided in Table 4. 
The estimated unit costs are US$ 2.695/m’. 

In Table 5, the unit costs are represented as 
percentages of the total unit cost, where capital costs 
account for approximately 65% of the total cost and 
operating cost is estimated to be around 35%. 

Table 5 shows the estimated unit costs using a small 
experimental plant operating data. A scale-up 
approach will be used to determine the cost of a large 
capacity plant. Due to economies of scale, the unit 
cost for a larger plant is likely to be lower than in a 
smaller one. The Eq. (1) will be used to scale-up the 
unit cost of the pilot plant: 

Cx = Cy Q/Q))" (1) 
where: 

C, = The capital cost for a large plant with a 
specific capacity; 

C, = The capital cost for a small pilot plant with its 


actual capacity; 
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Fig.1 Integrated submerged membrane microfiltration treatment system. 
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Table 1 Operating data for the integration system and the experimental plant. 


Parameter 

Plant capacity (m*/day) 
Plant life (years) 
Running days 

Running hours (h) 
Feed inflow (l/h) 
Filtrate outflow (l/h) 
Feed pressure (bar) 
Total feed (m°) 

Total filtrate (m°) 
Total electricity consumed (kWh) 


Table 2 Estimated costs of capital and operating expenses. 


Parameters 

Total unit capital cost 
Machinery—compressor 
Machinery—feed pump 
Machinery—transfer pump 
Machinery—backwash pump 
Tanks 

PVC pipes and valves 
Membranes 

Civil work 

Concrete base 

Kirby shade 

Operating cost 

Manpower for monitoring (2 technicians) 
Chemicals agent (100 Its) 
Maintenance and spare parts 


Dairy operations 
24 

20 
137 
2,972 
1,000 
1,000 
1.05 
2:972 
2972 
27.8 


Cost (US $) 
95,000 
3,000 
1,000 
900 
1,000 
2,830 
540 
9,260 
3,330 
3,330 
2,660 
0.335 
2,330 
400 
1,066 


Table 3 Total costs of dairy processing wastewater effluent for reuse. 


Capital cost 
Machinery—unit cost 


2 Machinery—compressor 

3: Machinery—feed pump 

4 Machinery—transfer pump 

3 Machinery—backwash pump 

6 Membranes 

7 Tanks (3 tanks) 

8 Manpower 

6 Total civil work 
Total capital cost 

B. Operating cost 

l Electricity = (E/Q)*R 

2 Chemicals 

3 Pipes and valves 

4 Maintenance and Maintenance parts 
Total operating cost 

Total cost 


Cost (US $) 
95,000 
12,000 
4,000 
3,600 
4,000 
37,335 
2,830 
4,665 
6,000 
169,165 


1.65 
1,200 
545 
1,065 
2,810 
171,975 


Cost Evaluation of Compact Dairy Wastewater Treatment System in Kuwait 115 


Table 4 Estimates of unit costs for the treatment of wastewater produced from dairy production in a small experimental 


pilot plant. 


Wastewater treatment 


A Total unit capital cost (depreciation) 
l Unit Cost 

2 Machinery—compressor 

3 Machinery—feed pump 

4 Machinery—transfer pump 

J Machinery—backwash pump 

6 Membranes 

7 Tanks 

8 Manpower 

9 Total civil work 

B Operation cost 

l Electricity 

2 Chemicals 

3 Pipes and valves 

4 Maintenance and maintenance parts 
Total cost (A+ B) 


Cost (US $) 


0.543 
0.068 
0.023 
0.020 
0.023 
0.212 
0.016 
0.811 
0.034 


0.0006 
0.403 
0.182 
0.359 
2.695 


Table 5 Distribution of unit costs for the treatment of wastewater produced from dairy production in a small experimental 


pilot plant. 


Wastewater treatment 


A Total unit capital cost (depreciation) 
l Unit cost 

2 Machinery—compressor 

3 Machinery—feed pump 

4 Machinery—transfer pump 

5 Machinery—backwash pump 

6 Membranes 

7 Tanks 

8 Manpower 

9 Total civil work 

B Operation cost 

1 Electricity 

2 Chemicals 

3 Pipes and valves 

4 Maintenance and maintenance parts 
Total cost (A+ B) 


Cost (%) 
64.90% 
20.12% 
2.54% 
0.85% 
0.76% 
0.85% 
7.85% 
0.60% 
30.06% 
1.27% 
35.10% 
0.02% 
14.98% 
6.78% 
13.32% 
100% 


Table 6 Unit cost estimates for treatment of wastewater produced from dairy production in a large commercial size plant. 


Unit cost (US $) 


n = 0.95 
Capital 0.63 
Operation 0.945 
Total 1.575 


Q, = The capacity of a large plant (120,000 m’°/d); 
Qy = 


n = A parameter representing economies of scale. 


The capacity of a small pilot; 


Economies of scale 
n = 0.90 


0.413 
0.945 
1.355 


n = 0.85 
0.270 
0.945 
L575 


The value of n is unknown due to lack of relevant 


information. Therefore, different assigned values will 
be used in this study (n = 0.95, n = 0.90 and n = 0.85) 
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which imply a modest to reasonable level of 
economies of scale. Not all cost components are 
affected to the same degree by plant capacity. The 
most affected component is the capital cost, so for 
simplicity, an assumption will be made that no 
economies of scale exist in other components. 

Table 6 shows the variation of the unit cost with 
different value of n .Whereas the economies of scale 
increase, the estimated unit cost decreases. 

For a modest economies of scale (n = 0.95), the unit 
cost of the system is estimated to be US$ 1.5/m? and 
for a greater economies of scale (n = 0.85), the unit 
cost is estimated to be US$ 1.575/m’. 

Currently, 17 tankers with an 30 m° capacity are 
hired to carry the wastewater. The cost of each tanker 
is US$ 45 per day. So the annual cost for wastewater 
discharge is a total of US$ 280,000. The total amount 
of wastewater to be carried for one year is 175,200 m’? 
and therefore, the annual cost of transporting 
generated wastewater is US$ 1.55/m°. Extra cost 
needs to be added if the tankers dispose the 
wastewater into nearby conventional treatment plant. 
The cost of conventional treatment of municipal 
wastewater in Kuwait was reported by previous study 
to be US$ 0.55/m° [7]. Therefore, the total cost of this 
process will be US$ 2.1/m’. 


5. Conclusion 


The unit costs of dairy-processing wastewater 
treatment for a small experimental plant using the 
integrated system technique have been estimated 
based on the operational data and cost elements. The 
result is US$ 2.695/m? which, in capital costs account 
for 65% of the total cost, whereas the operating cost 
accounts for 35% of the total. Using economies of 
scale to calculate the unit cost for a large commercial 


plant, the analysis shows that the unit cost is 


US$ 1.575/m° , whereas the unit cost for wastewater 
transportation and treatment in a conventional plant is 
US$ 2.1/m°. This reveals that the proposed treatment 
system is cost effective and can be used to treat 


industrial waste effluent. 
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Abstract: Genetic structure data of five populations of the Luehea divaricata Mart. & Zucc., forest tree species under development 
in the Atlantic Forest biome, obtained by microsatellite DNA markers, were used in simulations to study their reproductive and 
ecological pattern. Different selfing and migration rates were tested, using the observed and expected heterozygosity of 0.55 and 0.67, 
respectively, obtained through the use of microsatellite markers. Closest values were obtained with the use of selfing rates of 0.3 and 
migration of 0.2. These results suggest the presence of some self-incompatibility system between these species, which reduces, but 
does not prevent the self-fertilization. The migration rate contributes to a low genetic differentiation between the populations, making 
the reproductive mode, responsible for the inbreeding observed in the same populations. Authors suggest continuous monitoring of 
the genetic variability as a guarantee for the persistence of these populations. The study focus on the importance of using computer 
simulations to investigate ecologic, reproductive and genetic patterns for forestry populations, thus enabling the application of 
suitable measures for conservation. 


Key words: Computational simulations, conservation biology, inbreeding, heterozygosity. 


1. Introduction fragments, the reproductive system and distance of 


, , dispersion of pollen and seeds, the geographic 
The exponential growth of the human population, , , , o , 
, isolation can result in the reproductive isolation, and 
which has been observed over the last decades, has a , , a 
, , consequently restriction in the crossings inside the 
been accompanied by the destruction of the planet’s Eaa 
, , , fragments. In a long term, speciation may occur as 
diverse ecosystems, including forestry. The l 
S members of separate groups can develop different 
overexploitation of forest resources and the , , , 
l , evolutionary mechanisms to the point of no longer 
devastation of forest areas to the expansion of , , , 
, E , being able to cross, even if the gene flow is restored. 
agricultural activities are among the main human a ee 
2. However, the extinction of the population is most 
activities that have resulted in the destruction and 
D likely, due to the size of the isolates, since in small 
fragmentation of habitats and hence the extinction of , , , 
, populations there is a greater performance of genetic 
forest species [1]. , , 
, , : drift and a greater occurrence of inbreeding. These 
Habitat fragmentation is a phenomenon that , ae 
S ; events reduce genetic variability and make the 
consists of the subdivision of an area, geographically , , , 
, , | population susceptible to the effects of cyclical 
isolating members of a certain population into two or 
l , environmental changes such as the occurrence of pests 
more fragments. Depending on the distance between , 
or diseases, occurrence of frosts or floods [1, 2]. 
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evolution acts on an existing variability in the 
population, allowing the natural selection to preserve 
the individuals better adapted to a wide range of 
environmental variability [3]. 

The management of a fragmented population goes 
through the conservation of its genetic variability that 
can be accomplished by expanding the areas of its 
adds 


variability or by allowing the genetic flow between the 


fragments with germplasm that genetic 
fragments. In order to define a priority population for 
conservation and the conservation method to be 
adopted, the parameters of the genetic structure must 
be determined through the use of molecular markers, 
from which the rate of crosses and the level of genetic 
flow among its fragments can be estimated [4]. 
Computational simulations offer more accurate 
alternatives for the estimation of population genetic 
parameters, including crossover and gene flow rates. 
The EASYPOP program [5] was used in the present 


study to simulate selfing and migration rates from 











microsatellite markers data obtained from five 
populations of Luehea divaricata Mart. & Zucc. 
(Malvaceae) forest tree species, under development in 
a Brazilian area of the Atlantic Forest biome, in order 
to provide credible information for planning their 


conservation. 


2. Material and Methods 


2.1 Research Area 


Authors studied five Luehea divaricata Mart. & 
Zucc. (Malvaceae) forest populations located in the 


Brazilin States of Santa Catarina, Paraná, São Paulo 
and Minas Gerais (Fig. 1). 


2.2 Model Settings 


Authors used different selfing (0.1, 0.3 and 0.5) and 
migration rates (from 0.1 to 0.9, with steps of 0.1) in 
simulations with the computer program EASYPOP 


version 2.0.1 [5] to determine those that would result 
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Fig. 1 Geographic location of the five L. divaricata populations used in molecular characterization and simulations [6]. 
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Table 1 Values of observed and expected heterozygosities from the five populations used in the study. 


Population State A N 
Lages (LG) SC 2 32 
Irati (IR) PR 1.3 30 
São Jerônimo (SJ) PR 2.5 30 
Taciba (TA) SP 6.2 28 
Lavras (LA) MG T3 30 
Average 


Ho He Fis 
0.64 0.64 -0.01 
0.65 0.67 0.03 
0.54 0.63 0.15 
0.49 0.7 0.31 
0.43 0.71 0.40 
0.55 0.67 0.88 


A: area size, in hectares; N: number of individuals; Ho: observed heterozygosity, obtained by the use of microsatellite markers; He: 


expected heterozygosity, obtained by the use of microsatellite markers; Fig: inbreeding coefficient. Adapted from Ref. [6]. 


in parameters similar to those determined by Conson 
[6] when studying the genetic structure of nine Luehea 
divaricata populations in the Atlantic Forest, in Brazil. 
In the study, the 


heterozygosity was 0.55, and the average expected 


referred average observed 
heterozygosity was 0.67 (Table 1). 

For the simulations, authors considered diploid 
hermaphrodite specie, with non-random mating and 
no clonal reproduction. A spatial migration model was 
considered, with coordinates based on the geographic 
location of the populations and where obtained from 
the real geographical coordinates (latitude and 
longitude) according to the map. 

With respect to the mutation settings, authors 
assumed 10 loci evolving according to the single-step 
mutation model (SSM), with a proportion of 0.1 
K-allele model (KAM) events, under 50 possible 
allelic states, according to the data obtained by 
Conson [6]. A mutation rate of 0.0001 mutations per 
locus per generation was assumed. The genetic 
variability of the initial population was considered the 
maximum, and authors simulated 400 generations. For 
each combination of selfing and migration rates, 
authors replicated 100 replicates. 

The observed and expected heterozygosities (0.55 
and 0.67, respectively), from Conson [6] were used to 
select the model settings that presented values closest 
to the field observations (Table 1). Authors used the 
two independent sample t-tests at a 5% level of 
probability to compare the values of observed and 
expected heterozygosities of the selected model with 


others. 


2.3 FSTAT Analyses 


Authors used FSTAT program [7] to analyze the 
output files of the selected model. Through the use of 
FSTAT, the forest populations, including allelic 
richness, genetic diversity, and indices of genetic 
differentiation among populations were characterized. 
Authors estimated gene flow by the average number 
of migrants per generation, based on the proportion of 
the total genetic variance in a subpopulation (Rst). Rst 
statistic was used in this study in place of the common 
Fsr statistic because authors used microsatellite loci 
data that evolve according to step-by-step mutations 
and is the most indicative estimate for genetic 
differentiation between populations. 

The coefficient of inbreeding within the populations 
(Fis) was used to estimate the rate of apparent 
outcrossing (fa), given by the equation of Nei and 
Syakudo [8]: 

1 — Fis 
‘a= 74 Fis 


In turn, authors estimated gene flow in terms of the 


(1) 


average number of migrants per generation (Vm) from 
the value of Rsr [9]: 
1 


nm =z- (2) 


3. Results and Discussion 
3.1 Selection of Model Parameters 


Values of observed (Ho) and expected (He) 
EASY POP 


simulations for different selfing and migration rates are 


heterozygosities generated from 
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Table 2 Estimates of observed (Ho) and expected (He) heterozygosities generated by simulations, for different rates of 


selfing and migration 


Migration rate 0.1 
Ho He Ho 

0.1 0.656 0.694 0.543 
0.2 0.668 0.706 0.555 
0.3 0.670 0.708 0.565 
0.4 0.674 0.711 0.568 
0.5 0.681 0.719 0.571 
0.6 0.683 0.719 0.573 
0.7 0.684 0.722 0.575 
0.8 0.687 0.723 0.578 
0.9 0.689 0.725 0.580 


Rate of selfing 
0.3 0.5 

He Ho He 

0.658 0.424 0.633 
0.674 0.433 0.647 
0.684 0.433 0.649 
0.686 0.437 0.653 
0.692 0.441 0.659 
0.695 0.442 0.663 
0.698 0.444 0.664 
0.701 0.446 0.670 
0.705 0.452 0.674 


presented in Table 2. There is clear trend of increasing 
heterozygosity (both Ho and He) with increasing 
migration rate and reducing selfing. Heterozygosity is 
one of the parameters of genetic variability. According 
to Madsen et al. [10], gene flow enables the allele 
exchange between populations, thus increasing the 
genetic variability. In turn, selfing limits the pollen 
and seed dispersion and the potential for 


recombination between alleles from different 
individuals; as a result, autogam species tend to 
present low variability within populations and high 
variability between populations [11]. 

The closest values of Ho and He to those obtained 
through the use of microsatellite markers (0.55 and 
0.67, respectively) were generated under the selfing of 
0.3 and the migration rate of 0.2 and these values were 
significantly different from the other rates, based on 
two independent sample t-tests at a 5% level of 
probability of error. This model enables to 


characterize the populations. 
3.2 Genetic Structure of the Forest Populations 


The rate of selfing of 0.3 generated from 
simulations implies an outcrossing rate of 0.7. This 
value is in accordance with the apparent outcrossing 
(ta) of 0.696 obtained by using Eq. (2), given the 
intra-population inbreeding coefficient (Fis) of 0.179 
generated from the FSTAT analysis. This finding 


enables us to classify the mode of reproduction of the 


population as mixed, with a predominance of 
outcrossing, according to the classification of Destro 
and Montalván [12]. This author classifies populations 
as autogamous (outcross from 0 to 0.05), mixed 
(outcross from 0.05 to 0.95) and allogamous (outcross 
from 0.95 to 1.0). 

Mixed reproduction system in hermaphrodite 
vegetal species, like Luehea divaricata, suggests the 
presence of self-incompatibility mechanism in vegetal 
species that reduce, but do not completely prevent 
selfing [13]. 


self-incompatibility is frequently used by angiosperms 


According to Bawa [14], 


to prevent the adverse effects of inbreeding and the 
subsequent loss of genetic variability. Studying Ceiba 
aesculifolia populations, Quesada, et al. [15] observed 
outcrossing rates close to 1.0 because of a 
self-incompatibility mechanism present in this 
species. 

With respect to the migration settings, the 
simulations derived a rate of 0.2 between populations, 
which is enough to high genetic 
differentiation between (Rsr = 0.0062, from the 
FSTAT analyze). For gene flow, in terms of Nm, the 


prevent 


number of migrants per generation, determined 
through the use of Eq. (3), was 40. This is a high value 
according to the scale set by Govindaraju [16], which 
distinguishes the following levels for gene flow: high 
(Nm > 1), intermediate (0.25 < Nm < 0.99), and low 
(Nm < 0.25). 
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Fig. 2 Simulation of observed (Ho) and expected (Hs) 
heterozygosities across 400 generations in five populations 
of Luehea divaricata. 


Gene flow is the movement of genes (gametes, 
propagules, and individuals) in populations that 
effectively changes the spatial distributions thereof. 
High rates of gene flow result in low genetic 
differentiation between populations by opposing the 
effects of evolutionary forces (inbreeding and genetic 
drift) [17]. A study regarding genetic structure for 
Luehea divaricata in Brazil [9] showed high genetic 
differentiation (Fst = 0.22) between populations, due 
to low gene flow. According to Neigel [17], gene flow 
is an evolutionary force opposite to genetic drift, 
acting toward the homogenizing allele frequencies 
between populations, thus resulting in low genetic 
differentiation between them. 

The behavior of 


heterozygosity were plotted across the 400 generations 


observed and expected 
simulated (Fig. 2). High levels of heterozygosities 
were observed (values close to one) in the beginning, 
with a progressive loss generation by generation. 
Diversity loss may occur due to inbreeding, despite 
the self-incompatibility mechanisms present in this 
species. As a result of genetic drift effects, some 
alleles may lose every generation, including those for 
self-incompatibility. Inbreeding can also occur from 
biparental crossings when population size is small. 
According to Ridley [18], while a great habitat area 
may have sustained a single great population, it is 
possible that none of its fragments is able to sustain 


enough a subpopulation for a long time. 


3.3 Implications for Conservation 


Forest habitats have been destroyed and fragmented 
due to anthropic activities that turn a great and 
continuous area into two or more small sized 
fragments. In this condition, ecologic and genetic 
and lead 


consequences [1, 2]. Forest fragmentation affects the 


processes may occur to disastrous 
phenology, pollination patterns and reproductive 


success of species [2]. The reduction on the 


population size may reduce the density of 
reproductive trees and limit pollen availability, and 
propitiate the occurrence of inbreeding and genetic 
drift effects that lead to progressive loss of genetic 
variability and, ultimately, the population extinction. 
The findings in this study suggest five genetically 
stable populations with moderate losses of variability 
due to high rates of gene flow between them despite 
high selfing rates. Gene flow enables sharing alleles 
from different populations and increases the 
heterozygosity. In order to maintain the present 
scenario in the future, authors suggest keeping 
monitoring the genetic variability as a guarantee for 


the persistence of these populations. 


4. Conclusion 


The five populations of Luehea divaricata present 
moderate levels of inbreeding as a result of the 
reproductive mode in this species. However, high 
rates of gene flow between populations minimize the 
losses of genetic variability. 

This study stressed the importance of using 
investigate ecologic, 


computer simulations to 


reproductive and genetic patterns for forestry 
populations, thus enabling the application of suitable 


measures for conservation. 
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